Fast Numerical simulations of 2D turbulence using 
a dynamic model for Subgrid Motions 

J-P. Laval 1 , B. Dubrulle 2 ' 3 and S. V. Nazarenko 4 
February 2, 2008 

1 CEA/DAPNIA/SAp, CE Saclay, 91191 Gif sur Yvette Cedex, France 



NCAR, P.O. Box 3000, Boulder CO 80307-3000 

CNRS, URA 285, Observatoire Midi-Pyrenees, 14 avenue Belin, F- 
31400 Toulouse, France 



4 Mathematics Institute University of Warwick COVENTRY CV4 7AL, 
UK 

Abstract 

We present numerical simulation of 2D turbulent flow using a new 
model for the subgrid scales which are computed using a dynamic equa- 
tion linking the subgrid scales with the resolved velocity. This equation 
is not postulated, but derived from the constitutive equations under 
the assumption that the non-linear interactions of subgrid scales be- 
tween themselves are equivalent to a turbulent viscosity. This results 
in a linear stochastic equation for the subgrid scales, which can be 
numerically solved by a decomposition of the subgrid scales into local- 
ized wave-packets. These wave-packets are transported by the resolved 
scale velocity and have wavenumbers and amplitude which evolve ac- 
cording to the resolved strain and the stochastic forcing. The perfor- 
mances of our model are compared with Direct Numerical Simulations 
of decaying and forced turbulence. For a same resolution, numerical 
simulations using our model allow for a significant reduction of the 
computational time (of the order of 100 in the case we consider), and 
allow the achievement of significantly larger Reynolds number than the 
direct method. 
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1 Introduction 



The dynamics of high Reynolds number turbulent flow couples a large range 
of scales from the characteristic size of the domain to the dissipative scales. 
This range is usually too large to be fully resolved by Direct Numerical 
Simulation of the Navier-Stokes equations (DNS). This is the case for many 
applications like in aeronautics, geophysics or astrophysics where the typical 
Reynolds numbers are of the order of 10 6 to 10 12 . The largest Reynolds 
numbers which can be achieved by DNS are of the order of 10 4 to 10 5 111, 
0]. The difficulty of direct simulation of turbulent flows at high Reynolds 
number arises because of two scaling laws: first, the memory requirement to 
be able to deal with all the scale of motions growths like -Re 9 / 4 , where Re is 
the Reynolds number; second, the time stepping to be used to advance the 
equations has to be computed using the smallest resolved scale. As a result, 
computational time usually growths like Re 5 . Despite of the growth of the 
computer power ||, the direct simulation of all realistic problems will not 
be conceivable before a long time. 

Part of these difficulties could be circumvented if one could find an ef- 
ficient approximate scheme to model the small scales of the turbulent flow, 
which usually monopolize 90 per cent of the computations and, for most 
applications, do not need to be known with the same precision as the large 
scales. From an algorithmic point of view, several methods have been pro- 
posed to try to describe the small turbulent scales with less degree of freedom 
than in a DNS like for example, sparse Fourier transform pL rarefied modes 
||, wavelets Q, self-adaptative grid mesh 0, |8). From a theoretical point 
of view, the effort has been put on the model of the action of the small 
scales onto the larger scales of motion, or onto the small scales themselves, 
to be able to deal with simpler and cheaper systems to compute. In the ex- 
treme case of large eddy simulations (LES) or random averaged numerical 
simulations (RANS), the effect of the small scales at large scales is directly 
modeled as a function of the resolved or mean scales of motions, resulting 
in both memory and computational gain. This is often done at the price of 
introduction of arbitrary parameters, which require calibration and which 
lower the predictive power of the simulation. Also, some models are not 
theoretically satisfying, since they break for example the basic symmetries 
of the original equations || . 

In this paper, we continue a study initiated in [ 10 1 and develop a new 
dynamic model for the turbulent subgrid quantities which is directly de- 
rived from the constitutive equations of motions. This model respects all 
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the basic symmetries of the original equations, as well as their conservation 
properties. The main assumption underlying the model is that the smaller 
turbulent scales are much less energetic than the resolved scales, and so, 
that their main interactions are with the resolved flow. It then appears rea- 
sonable to try to model these main interactions with as few approximations 
as possible, while more freedom is allowed regarding their mutual interac- 
tions. In the model we consider, the interactions between the resolved and 
subgrid scales are taken into account without any approximation, while the 
mutual interactions between subgrid scales are replaced by a turbulent vis- 
cosity. In 2D geometry, where energy condensation at large scale guarantees 
that the dynamics of the small scales is mainly non-local in scales this tur- 
bulent viscosity is so small that it can be set to zero in some simple cases 
[11|. This allows a parameter free model of the subgrid scales, which will 
be described in the first part of this paper. In situations where dynamic 
processes at small scales, like vortex stretching, are important to regulate 
any small-scale instability, or even in 2D when "large" and "small" scales 
are not well enough separated, the turbulent viscosity cannot be set to zero 
[12|. This option introduces a free parameter into the model, but, since it 
appears at a sub-dominant level, we may expect that its choice is not so 
critical to the success of the model. 

After parameterization of the non-linear interactions between the sub- 
grid scales, their resulting equation of motion becomes linear. This essential 
feature has two advantages: first, it enables, in certain simple flow geome- 
tries, analytical solutions for the subgrid-scale dynamics as a function of the 
resolved quantities, thereby closing the equations of motions at the resolved- 
scale level. This property was used to obtain analytical solution for mean 
profiles in channel flows |l3|, [l4|, |l5[ or for the Planetary Surface Layer 
]j~6| , |i~7| |. Second advantage of the linear description is that it allows for 
more efficient algorithms of integrations, using Lagrangian methods where 
the time stepping is done via a criterion based on the resolved scales. This 
allows numerical computations with a larger time step than DNS, and thus, 
at a lower computational cost. 

An essential part of our model is the averaged Reynolds stresses de- 
scribing the feedback of the subgrid scales onto the resolved component. In 
[fLOf , we considered the scale separation parameter to be much less than the 
nonlinearity of the small scales and derived a feedback term in which the 
quadratic in the small-scale amplitude terms gave the main contribution into 
the Reynolds stress. Then, the model that consisted of two coupled "flu- 
ids" (the resolved and the subgrid one) was used to solve several problems 
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with ([jO]]) and without scale separation ([[Top), the later being the typical 
problems such as the forced turbulence, the vortex merger and a turbulence 
decay. It was noticed that in the problems without a natural scale sepa- 
ration, such a model described very well the small-scale dynamics whereas 
improvements in modeling of the large scales (which is the main aim of the 
LES) was much more modest. The situation was clarified in [11] where 
direct numerical evaluations of the different contributions to the Reynolds 
stress were made. It was shown that in the problems without the scale sep- 
aration the dominant contribution to the Reynolds stress comes from the 
linear rather than quadratic in the small-scale amplitude term. 

In the present contribution, we focus on advancing the numerical ap- 
proach introduced in [ 10 1 by introducing a better model for the turbulent 
Reynolds stress based on the results of the a priori tests performed in [11]. 
Also, a more efficient procedure is used in this paper to generate the subgrid 
vorticity and velocity fields. Our goal here is to test our model both as a 
numerical method for an improved DNS, in the sense that we can compute 
the whole range of scales at a lower computational cost, and as pseudo LES 
method, in the sense that we compute only the larger scales, at an even 
lower computational cost. For the sake of simplicity, we shall consider here 
only the two-dimensional case, where both the hypothesis of the turbulent 
model [11] and its consequences [^, 10] have been studied in detail. The 
generalization to the 3D case is the subject of an ongoing research. 



2 The turbulence model 

2.1 Scale decomposition 

We consider a two-dimensional incompressible inviscid flow obeying the 
equations: 

dtu) + div(uw) = vAu, , , 

divu =0, ^ ' 

where u is the velocity and v is the viscosity. In 2D geometry, the vorticity 
has only one non-zero component that we denote by oj. The resolved- and 
the subgrid-scale part of the velocity and the vorticity fields are defined via 
a filtering procedure: 

U(x,t) = u(x,t)= f G(x-x')u(x',t)fix / (2) 
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n(x,t) = tj(x,t) =J G(xL-x r )u(x r ,t)dx.' (3) 

Here, G is a filter, such that the resolved scales contain the main part of the 
total energy. In Section III.C, we shall propose a special shape of G obeying 
this condition, and which is well adapted to our numerical method. 
Each field is then decomposed as follows 

u(x,t) = U(x,t) + u'(x,t), (4) 
w(x,i) = n(x,t) +u/(x,i), (5) 

where the primes denote subgrid-scale quantities. Inserting this decompo- 
sition into the Navier-Stokes equations and separating the resolved and the 
subgrid-scale parts, we get a set of coupled equations: 



d t fl + div (Ufi) + div (u'fi) 

+ div (UuT) + div (uV) = ^A$7, (6) 



and 



«W + div (UO) - div (Ufi) 

+ div (U«/) - div (TO) 

+ div (u'n) - div 

+ div (uV) - div (t?ci7) = i^Aw'. (7) 

The second up to the sixth terms of the l.h.s. of ([?]) are the contributions 
due to non-local interactions between the resolved and the subgrid scales. 
The last two terms in the l.h.s. of (|^) are the contributions due to non- 
linear (local in scale space) interactions among the subgrid scales. In the 
favorable case where most of the energy is at the resolved scales, these last 
two terms can be expected to negligible with respect to, e.g., the second 
up to the fourth terms of the l.h.s. of (0). Indeed, in a recent detailed 
numerical analysis of the system (R3) and (|?]) , Laval et al showed that for 
a very steep filter (when the filter G is a cut-off in the spectral space), the 
small-scale dynamics (and thus the large-scale dynamics) is independent of 
the local interactions: when the later are neglected, the small-scale velocity 
and vorticity field are not significantly changed, in both forced and decaying 



turbulence, even after several eddy turn-over times [11]. The analysis has 
also shown that in 2D turbulence, the leading order contribution at both 
large and small-scale comes from the correlations involving the large-scale 
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velocity field, i.e. and Uw'. The next leading order contribution comes 
from the correlations between small-scale velocity and vorticity u'uj', while 
the correlation between small-scale velocity and large-scale vorticity u'Q 
gives the lowest contribution. 

2.2 A priori numerical estimates 

In the present case, the filter G is smoother, and subgrid scales include both 
large and small scales. It is thus important to conduct additional numerical 
evaluation of the various terms of (|6|) and ([?]) using a smooth filter, to 
check the validity of the non-locality assumption for the subgrid scales. We 
performed an apriori test with the filter which will be used in our model. 
The choice of the filter will be discussed in section [D|. The test was done 
using a vorticity field from a DNS on a 1024 2 grid. The Figures |l] and || 
give the distribution of each non-linear terms involved in the resolved scales 
and the subfilter-scales equations. The result of the filtering process one the 
initial field is shown fig. |3|. Even with a smooth filter, the non-linear term 
involving only subgrid scales are still small compared to the higher order 
term. 

2.3 The model 

Keeping only the leading order contributions, the coupled equations @ and 
([?]) become: 

d t ft + divTjn + divUZ7 = i/Afi, (8) 



d t J + div (Uu/) = F + v t Au/, 

F(x,t) = -(div(UO) -divUH) + divlI7. (9) 

Here, vt is a turbulent viscosity, which will be introduced to damp the small- 
scale noise arising in our numerical scheme and F is a force which describes 
the subgrid-scale generation via the enstrophy cascade (energy cascade in 
3D). 

3 Numerical implementation 
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Figure 1: Comparison of the square moduli of each component of the non- 
linear part of the large scale equation || in Fourier space. The filter used for 
the scale separation is a smooth filter in Fourier space defined by eq. fl23|). 
The energy spectra of both resolved scales and subfilter scales are shown fig. 

I 

3.1 Numerical strategy and performances 

The resolved-scale equation has the form of an Euler equation with an ad- 
ditional forcing coming from the interaction with the subgrid-scale motions. 
One is therefore led to standard strategies to solve this equation, depending 
on the type of the flow (spectral methods for flows in simple periodic geom- 
etry, or finite difference or finite elements for more complicated geometries) . 
Here, we shall consider the periodic case, thereby using a spectral code for 
solving the resolved-scale equation. This code is described in |n| . 

The situation is markedly different at the subgrid level, where the ba- 
sic equation is linear in the subgrid motions, with an inhomogeneous part 
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Figure 2: Same comparison as in fig. [l] but now for the non-linear terms in 
the subfilter scale equation. 

provided by the subgrid scale generation via the enstrophy cascade. This lin- 
earity suggests a solution strategy based on projection of the subgrid scales 
onto appropriate modes. Since the subgrid scales are usually very inhomo- 
geneous (e.g., vorticity filaments), it seems logical to use a decomposition of 
the subgrid-scale velocity field into localized modes, thereby optimizing the 
memory requirements to store the subgrid-scale fields. A popular local mode 
decomposition uses wavelets (see e.g. |2(J). Wavelets are however sometimes 
difficult to implement, and are not very handy to use in analytical compu- 
tations. Here, we choose to use a Gabor decomposition, which provides a 
localized description while allowing theoretical manipulations similar to that 
obtained with Fourier modes. The linearity of the subgrid motions also pro- 
vides room for a further computational time reduction via semi-Lagrangian 
methods of integration, using a time step related to the resolved scale. The 




div(UQ)(k) - <div(UQ)>(k) | 2 
div(u'n)(k) - <div(u'Q)>(k) | ; 
div(Uw')( k ) - <div(Uc/)>(k) 
div(uV)(k) - <div(uV)>(k) 
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Figure 3: Energy spectra of the resolved and the subfilter components for 
the same run as in fig. [I] and ^. The filter used for the scale separation 
is defined in eq. (|23l ) with dh = 27r/32. The initial field is the result of a 
Direct Numerical Simulation of decaying turbulence on a 1024 2 grid. 



numerical expected performances of our model, integrated using the semi- 
Lagrangian, Gabor method, at a given resolution and total integration time, 
is given in Table 0, and compared with both standard spectral methods 
based on FFT, and a LES approach based on FFT. Computational time for 
several particular examples is given in table [2] which shows an obvious gain 
in speed for our method with respect to the DNS. If we now compare with 
a classical LES, our method is somewhat slower because it is necessary to 
keep a sufficient number of localized modes for accurate description of the 
Reynolds stresses. However, a substantial gain in accuracy is achieved in 
comparison with traditional LES because of the better description of the 
nonlocal interaction of scales. 

3.2 Description of the algorithm 

Our algorithm of resolution of the turbulent model is based on five steps: 

1 Compute the force F (||) using the resolved and subgrid fields at time 

t; 

2 Project this force onto a set of Gabor modes; 



9 



Model (resolution) 


integration time 


DNS [FFT code] (N x N) 
LES [FFT code] (M x M) 
Our Model [FFT + Lagr.] 


d * NHog(N) 
c 2 * M 2 log(M) 
c 3 * M 2 log(M) +c A *N p 



Table 1: Comparison of expected performance between DNS, Classical LES 
models (like APVM or HDNS) based on FFT on a (M x M) grid and our 
model on a (M x M) grid for fully resolved scales and N p modes for approx- 
imated subfilter scale equation (see section [T^). 



Model 


r = 15 


r = 50 


Comp. time 


DNS (N = 1024) 






~ 10 days 


MO (M = 64, N p = 512 2 ) 


0.99403 


0.89794 


~ lh30 


Ml (M = 64, N p = 512 2 ) 


0.99450 


0.82523 


~ lh30 


M2 (M = 64, N p = 128 2 ) 


0.99409 


0.89735 


~ 8 mn 


APVM (M = 64 2 ) 


0.95882 


0.77429 


~ 2 mn 


HDNS (M = 64 2 ) 


0.93893 


0.53231 


~ 2 mn 



Table 2: Correlation between large-scale (k < 21) vorticity from DNS 
and vorticity field from : our model with a maximum of 512 2 particles (M0 
and Ml) and 128 2 particles (M2), the APVM model and an hyper- viscous 
simulation both on the resolved-scale grid iV = 64 2 . The same viscosity as 
in DNS was introduced in the resolved scale equation for the Ml simulation 
whereas no viscosity was introduced for the M0 and M2 simulations. The 
correlation coefficient are given at two different times (15 and 50 turnover 
times). The computations time with a "Sun Ultrasparc 3000 workstation" 
are given in the last column. 
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3 Advance the subgrid field to time t + dt in the Gabor space using this 
projection and the Lagrangian advection algorithm; 



4 Compute the resolved Reynolds stresses at time t+dt using a procedure 
similar to an inverse Gabor transform; 

5 Advance the resolved velocity field at time t + dt using this Reynolds 
stress. 

Step 1 and 5 involve standard procedures linked with our resolved scale 
code (in the present case a Fourier spectral code with an Adams-Bashforth 
time stepping where a full desaliasing was introduced by keeping only the 



first 2/3 smaller wavenumbers in each direction, see [19|). Step 2, 3 and 4 
involve new original procedures, based on interesting properties of the Gabor 
transform. For the sake of clarity , we present here only main results, leaving 
detailed computations in appendices. 

3.3 The continuous Gabor transform 

The Gabor transform (GT) is defined as 

u'(x,k,i) = |/(e*(x-x'))e ik ( x - x ')u'(x') ( ix', (10) 

where / is a rapidly decreasing function at infinity. Note that 1/e* has 
a meaning of the scale separation length. It has to be chosen to lie in 
between of the integral scale and the Kolmogorov scale, the later being far 
in the subgrid range in many applications. For LES purposes, 1/e* has to 
be close to (but not less than) the minimal resolved scale (grid scale dh in 
our case). Thus, to derive our model we need to assume the subfilter-scale 
wavenumber k to be greater than e* and use e*/k as a small parameter. 
Technically, however, is convinient to keep k fixed and perform expansions 
in small e*. 

In the special case where / = VG, where G is the filter (see eq. (0)), we 
have the following simple reconstruction formulae: 

u '( x >*) = m ^2 1 t ,n\ I u'(x,k,t)dk, (11) 



(2tt)* /(0) 

W(x, t) = J p(x, k, t)/(-k, t)] dk + 0(e*), (12) 
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where 3ft means the real part, and finally 

1 



uX-(x,t) 



(2n? 



^(x,k,t)u'(x,-k,t) dk + 0(e* 



(13) 



One may also note the connection between the Gabor velocity and the Gabor 
vorticity (similar to the Fourier quantities): 



u(k,t) 



k 2 



u(k,t)(e z x k)+0(e*), 



(14) 



where e z is a unit vector in z-direction. 

Finally, one may show that the GT of the subgrid-scale equation 
(see appendix [A] and [13, ffq] for more details) 



is 



(d t + XJ-V x - V x (U • k) • V fc ) w'(x, k, i) = F(x, t) - i^fcV(x, k, t). (15) 

3.4 The discretized optimum Gabor transform 

The continuous Gabor transform is not suitable for numerical implemen- 
tation. Moreover, because it involves both position and wavenumbers, it 
theoretically describes a given field via a much larger mode numbers than 
in traditional spectral methods (typically N 2D versus N D in D dimensions). 
However, as discussed previously, inhomogeneous fields may be represented 
with a high precision by a much smaller set of Gabor modes than the theo- 
retical number, a number even potentially smaller than N D (see e.g. |2l|). 
The algorithm we devised uses N D Gabor modes and, therefore, there is a 
room for further code optimization. 



3.4.1 The discretization 

The discretization of the subgrid field in the Gabor space is done via a 
Particle In Cell (PIC) method. This method is traditionally used in plasma 
physics, but was used previously in hydrodynamics by Nazarenko et aJ (2^] to 
study the interaction of sound wave-packets with turbulence. Details about 
the method can also be found in flpfl . In this method, the Gabor modes 
of the subgrid-scale vorticity field are replaced with discrete wave-packets 
(particles). For example, the vorticity field u/(x, k, t) is discretized with N p 
wave-packets a in the Gabor space: 

N p 

w'(x,k,t) = S x (x-xa(t)) S^k-M*)). (16) 

a=l 
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In (|T(j), a labels the wave-packet. Note that our representation allows several 
wave-packets with different wavenumbers to be located at the same position 
x, like in the Gabor transform. Here 5 X and are some interpolating 
functions (particle "size" ) . Since the computation of non- linear terms in the 
resolved-scale equation involves a second order spatial derivative, we are led 
to choose a linear interpolation in the x direction. In the k direction, a 
zero particle "thickness" turns out to be sufficient for our purpose. We thus 
adopt the following representation, 

S x (x) = S(x)S(y), (17) 
S k (k) = <%)%), (18) 

where x = (x,y), k = (p,q), 5 is the Dirac function and the function S(rj) 
is defined by: 

S(ri) = { ^ dh ~ ^D/M if M < dh i H9) 
1 otherwise, 

(20) 

Here, dh is a length scale governing the accuracy of the discretization. Its 
choice will be discussed later on. 

The PIC algorithm of reconstruction of the continuous fields from the 
discretized coefficients provides a natural filter G and the function / to be 
used in the Gabor transform through (see appendix |D[): 

G(x) = / 2 (x) = / 2 (0) x <? x (x), (21) 

with 

7(0) = , = (22) 



In the k space, the filtering operation is achieved by simply multiplication 
of the Fourier coefficients by the Fourier transform of G, which is 

A, \ 36 \ ( sin{pdh) \ ( sin{q dh) \\ 

q) = WW^dW 1 1 " ~Vdh~ ) I 1 " -ilh- ) / (23) 
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3.4.2 Reconstruction of the subgrid scale correlations 



With our discretization, the formulae of reconstruction of the subgrid corre- 
lations (|l^) can be obtained using (|16). Taking into account the condition 
that uj' is real, we have (see appendix C]for details): 



N B 



Uu/(x,t)~2U ^k; + r(k Q+ )]5 x (x-x a+ ), 



(24) 



Q+ = l 



where the sum is only over the wave-packets with positive wavenumber pj^. 
Similarly, one can use ( |45|) and (16) to re-write ( |i~3| ) as (63) (see appendix 

0) 



No 



u'v'(x,t) = 2 



a + = 
N v 



2\2 



,' 2 (x,t)=2 £ 



a+ = 
N„ 



2\2 



„/2, 



x,*) = 2 £ 



1 (.Pa+ + Qa+ 



2\2 



| ^"x (■"■ "a . , 



(25) 



3.4.3 Optimum discrete Gabor transform 

For a given field u/(x, t), associated with a grid of size 2ir/N, our "optimum" 
discrete Gabor transform only retains one wavenumber per position (i.e. N 2 
Gabor modes, for a N 2 initial discrete field). This wavenumber, and the 
amplitude of the corresponding Gabor mode are chosen as follows: we use 
the exact relation between the field and its iV 2 wave-packet components 



li/(x, t) 



/(0) 



x x Q j . 



(26) 



«4- = l 



An easy way to satisfy (26) exactly is to create one wave-packet at each grid 
point, so that N p = N 2 . Applying equation ( [26] ) at each grid point Xj we 
then obtain: 



1 . Since the vorticity w'(x, t) is real, each wave-packet at position x a and with 
wavenumber k a will have a partner wave-packet with same amplitude, at the same lo- 
cation and with an opposite wavenumber — k a . 
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u'(Xi,t) 



/(O) 



^k(t)]. 



(27) 



This equation fixes the real part and the physical coordinates of the 
wave-packet. To find its wavenumber, we use the velocity at the grid point 
u'(xj) to find two conditions: 



u'(xj,t) 

v'(Xi,t) 



/(0) (pI+qD 

2 +Pa t 

/(0) (Pl+ql) 



3k(<)], 
sk,(t)]. 



(28) 
(29) 



These conditions link the two components of the wavenumber with the imag- 
inary part of the amplitude of the wave-packet, which is still a free parame- 
ter at this stage. We then select the phase of the GT by requiring that the 
imaginary part of the wave packet equals its real part. Other choices could 
have been made, but this one turns out to simplify the computations. The 
characteristics of each wave-packets on can finally be summarized as: 



r »w 



Pi 



-u' f (x.i)/v' f (xi), 

v' f {Ki){l + (q a J PoLi ) )' 
Xi. 



(30) 



This procedure creates iV 2 wave-packets, from any vorticity field on a grid 
of size 2tt/N and can be used to initialize the subgrid scale vorticity field 
from any given initial condition or to GT the force F by transforming it into 
an equivalent vorticity field w^(x, t) = F(x,t)/dt. 



3.5 The Lagrangian scheme 

A Lagrangian interpretation of the subgrid scale equation ( ]T5| ) shows that 
its time integration is equivalent to the evolution of each wave-packet in 
the (x,k) space. Each wave-packet carries a complex amount of Gabor 
vorticity (cr a ) and is advected at the resolved-scale velocity U(x), while its 
wavenumber and amplitude evolve according to the local resolved strain. 
For the trajectory of the wave-packet and its amplitude we have 

x a = U(x a (t)), (31) 
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k Q = -V x (k Q (t) • U(x Q (i))) , 
& a = F(x a ,k a ,t) - v t k 2 a o a . 



(32) 
(33) 



By these equations, the wave-packets evolve continuously in the physical 
space in between the grid points associated with the finite, resolved-scale 
resolution. To find the resolved-scale quantities at the position of each par- 
ticle (i.e. possibly in between grid mesh points), an interpolation procedure 
is used via the functions S x and associated with the PIC method. 



3.6 Initial conditions 

All the simulations used the same initial random vorticity field with all the 
energy concentrated at very large scales ( the initial energy spectrum is 
given by E(k) = ke~( k ~ ko ' 1 with k a = 1). In practise, the noise associated 
with the PIC method is a function of the number of wave-packets used for 
computation. If the initial condition is such that the initial grid on which 
the vorticity is defined is too coarse (M too small), one can use a two step 
procedure to reach a reasonable number of wave-packets: in a first stage 
of the simulation, we create at each time step M 2 wave-packets with the 
vorticity (jj(x, t) = F(x,t)/dt created by the forcing F. During this stage, 
the wave-packets which were created at earlier times are moved into the 
phase space using the ray equations (31 J32]), but their amplitude is kept 



constant. This procedure is used until the total number of wave-packet 
reaches a desired number. From this point on, the wave-packet creation is 
shut down, and the wave-packets are evolved according to ([31, 32], 33). 



3.7 Noise reduction and effective turbulent viscosity 

For inviscid simulations (when the turbulent viscosity is taken equal to 0), 
the numerical procedure develops a noise at very subgrid scales, which is 
a function both of the integration time and the total number of particles. 
We have found that an efficient way to reduce this subgrid-scale noise could 
be achieved by periodically recreating a new set of N p wave-packets via 
first a rebuilding of the vorticity field in the physical space, on a grid of 
size y/Np using (|TT|) , followed by a re-creation of wave-packets using (|30|) . 
This procedure keeps the correlations associated with the subgrid-scale field 
unchanged, and therefore does not affect directly the resolved scale. In 
effect, this procedure acts as a filter for the subgrid scales which are smaller 
than the size of the reconstruction grid, and thus, can be seen as an effective 
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turbulent viscosity, which cannot be estimated a priori but which adjusts 
itself to the noise level. 

3.8 Test of the accuracy of the method 

To test the accuracy of our numerical method, we have performed a priori 
tests using data from a DNS at high resolution. The total vorticity field was 
divided into a subgrid-and a resolved-scale field using the same filter than in 
our model. The subgrid-scale field is further discretized into wave-packets 
according to eq. (|30|). A comparison was then made between the various 
components of the stresses in the DNS, and in our discretization scheme. 
The result of is shown fig. |], which shows the square modulus of each of 
the three Reynolds stress components with respect to their wavenumbers. 
The accuracy of the discretized reconstruction of each terms is remarkable 
for the large scales. At very small scales (k > 42), the discretization tends 
to produce a noise in the components of the Reynolds stress involving the 
subgrid-scale velocity. In our case, we do not consider these terms so that 
we do not have to bother about this noise. In situations where one, or 
two of these troublesome terms are considered, the noise can be removed 
via a filtering at small scale. In spectral simulations, this truncature is 
provided naturally via the procedure to remove aliasing, which filters out 
modes larger than 2k max /3 where k max is the maximum wavenumber of the 
resolved-scale field. The inset of figure || shows the comparison between the 
total real Reynolds stresses versus the total modeled Reynolds stress after 
desaliasing. The agreement is nearly perfect. 
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Figure 4: Moduli of the 3 non-linear terms in the resolved-scale equation 
@ computed directly by the spectral method (lines) compared with the 

same fields rebuilt by the PIC method (symbols). Here, + and are used 

for the field |div (\Jio')(k)\ 2 , A and — • — for |div (u'u')(k)\ 2 and O and 

for |div (u'$7)(/c)| 2 ). Corresponding comparison of the true and the 

PIC-modeled Reynolds Stresses r = div (uw) — div (UO) at each grid point 
Xij is shown in insert. 
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4 Numerical results 



Our method can be used to perform two kind of simulations, depending on 
whether one is interested in the small-scale behavior or not. In the first case, 
one needs to model the subgrid scales with a large number of wave-packets, 
so as to be able to reconstruct the small-scale field with a good accuracy. 
Typically, one needs about N p = N 2 wave-packets to be able to reconstruct 
faithfully details at scale 2ir/N and to produce a result which may be com- 
pared with a Direct Numerical Simulation at resolution N 2 . In the case 
when only large scales matter, one needs to keep the smallest number of 
wave-packet necessary to compute accurately the Reynolds stresses at the 
wavenumber cut-off. Since the resolved- scale field at the cut-off k c = M/2 
produces (by non-linear beating) information up to scale 2-k/M, we used 
a minimum of N p = M 2 wave-packets in our "Large Eddy Simulations". 
Finally, note that our method allows nearly inviscid computations, since it 
does not require the existence of a viscosity at large scale and since it uses a 
minimal "effective viscosity" which starts acting only at scales 2n/M. Very 
large Reynolds number can then be achieved via an adequate number of 
wave-packets. 

We present results illustrating these points in two classical situations, - 
decaying and forced turbulence. In each case, the performance of our model 
are discussed and compared with results from the DNS and other popular 
subgrid-scale parameterizations used in 2D turbulence. 

4.1 Decaying turbulence 

For the decaying case, the initial condition was chosen so that the energy 
is concentrated at very large scales. The reference DNS was performed at a 
resolution iV 2 = 1024 2 with a viscosity v = 1.8 10~ 4 , leading to a Reynolds 
number Re ~ 10 4 . The simulation was stopped after approximatively 50 
turnover times, by which time the initial condition has evolved into a robust 
dipole structure. The separation between resolved and subgrid scale used 
for our model is taken at k c = 21, corresponding to a computation over a 
grid 64 2 . 

Let us first compare the DNS to our model using a large number of 
wave-packets (N p = 512 2 ). We ran two different simulations: one in which 
the viscosity at resolved scale was set to v = 1.8 1CT 4 , like in the DNS 
("viscous" simulation); another one in which u was set to zero at resolved 
scale ("inviscid" simulation). The total simulation time of each simulation 
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is about lh30 on a Sun workstation, roughly 150 times less than the DNS. 

The total vorticity field after 50 turn over times, reconstructed by adding 
the contribution from the resolved and subgrid scales, is shown in Fig. ^ 
and compared with the vorticity field of the DNS. A global comparison for 
small and large scales can also be done via energy spectra. This is done in 
Fig. §. 




Figure 5: Total vorticity field of decaying turbulence after 50 turnover times 
as computed by DNS on a 1024 2 grid with viscosity v = 1.8 10~ 4 (left) 
compared with the same field computed by the model MO described in table 
2 (right). 

Clearly, the two pictures display good overall similarities at large scales, 
and marked differences at smaller scales. The largest structures are very sim- 
ilar and they are well localized even after 50 turnover times. The spectra in 
the two simulations overlap. The viscous simulation gives large-scale struc- 
tures which are in closer agreement with the (viscous) DNS (see Table 2). 
Together, these results confirm that at large scale, the dynamically impor- 
tant coupling term between resolved and subgrid scale is the term div(Uu)'), 
in agreement with the dynamical analysis of Laval et al [O], performed us- 
ing a cut-off filter. At smaller scales, our model seems to produce thinner 
filaments and smaller structures. This effect is clearly visible on the spectra 
of the two simulations: in the DNS, the k~ 5 inertial law starts to level off 
towards k = 60 (due to viscous effects), while in our model, the power law 
extends over a wider range of scales, up to approximately k = 100. To check 
whether this difference comes from our subgrid scale scheme, we performed 
a direct comparison of the smallest scales of the simulation k > 21. Fig. 
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Figure 6: The energy spectra obtained by the same simulation as in fig. 2 
compared with the spectra obtained by some other methods including our 
model with parameters MO from the table 2. 



I?] shows such a comparison at some earlier time (after about 15 turnover 
times), when the small scales have not yet been washed out by viscosity. At 
this time, differences are negligible, proving that our model also captures 
the dominant coupling mechanism at subgrid scales. This finding is also 
in agreement with the dynamical analysis of Laval et al 11]. The further 
differences arising over longer time scales can therefore be due to two effects: 
one is the error accumulation due to sub-dominant neglected terms in our 
model (like the terms involving coupling of the resolved vorticity SI with 
the subgrid velocity u'); the second is viscous effects. We believe that the 
first possibility is ruled out by the dynamical analysis of Laval et al [11], 
in which numerical simulations of our model equations were performed us- 
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Figure 7: Small-Scale vorticity field of decaying turbulence after 15 
turnovers. Result from a DNS on a 1024 2 grid with v = 1.8 1CT 4 (left) 
and from our model (MO). 

ing the same viscosity than in the direct numerical simulations 0. In that 
case, the differences at small scales appear to be negligible. We therefore 
interpret the differences between the two models as a viscous effect, and, 
actually, as an indication of the lower effective viscosity in our model. This 
would explain both the later bending of the spectra, and the finer structure 
of the filaments. 

The second series of comparison were performed using the minimum 
number of wave-packet, equal in the present case to N p = 128 2 . In such 
a case, the number of wave-packets does not allow accurate reconstruction 
for the scales less than 27r/128. The simulation is however much faster, and 
takes only 8 minutes of computational time. The total large-scale vorticity 
field after 50 turn over times is compared with the corresponding large-scale 
vorticity field of the DNS in Fig. |8[ The agreement is still very good and 
to quantify this agreement, we computed the correlation coefficient between 
the two simulations, defined by: 

n El< 4 J<64 w l( i >J>2(i,j) 

Y / Ei<i,i<64^i(^i) 2 Ei<ij<64^i(^i) 

where ui\(i,j) and uj2(i,j) are the large-scale vorticity field. This corre- 
lation coefficient was computed at two different times, corresponding to 15 

2 These simulations were not fast: they were based on spectral methods and were even 
more slowly than the direct method. 
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and 50 turnover time, and reported in table ^. At earlier time, our three 
models (viscous or not) are all characterized by a very good correlation co- 
efficient (about 99 percent). At later time, a slight difference appear. In 
fact, the best correlations are achieved by the viscous model, and the model 
with the least number of subgrid-scale modes. This is not surprising, since 
our noise removing procedure produces an effective viscosity which is larger 
as the number of subgrid-scale modes decreases. The high resolution model 
Ml is then the less viscous model of all three, and therefore, can be expected 
to produces the largest difference with respect to the viscous DNS. 

Two other popular 2D turbulent model were also tested along the same 
line. In the first one (HDNS), the viscous term vAQ of the Navier-Stokes 
equation is replaced by a "hyper- viscous" term v p A p £l. In our simulations, 
we took p = 8 and v p = 10 -18 . The second model is the Anticipated Poten- 
tial Vorticity Model (APVM) developed by Sadourny and Basdevant f23[| . 
Both models are very cheap, taking only about 2 minutes of computational 
time. They are also less accurate, as can be seen from both the Fig. ||, 
and the table. In the HDNS 64 2 , the correlation coefficient is only about 50 
percent in the end of the simulation. For the APVM, it is higher (about 75 
percent), but still lower than our "minimal model". The energy spectra of 
the 4 simulations can also be compared. This is done in Fig. ^. Our model 
develops an energy spectrum very close to the APVM one, and slightly less 
steep than the DNS at scales close to the cut-off. This is because at this 
scale, viscous effects start being felt and tend to bend the spectrum. The 
HDNS spectrum is much steeper near the cut-off than both the reference 
DNS and the two other models. 

4.2 Forced turbulence 

Similar simulations were performed in the case of forced turbulence. In such 
a case, our model was ran with M 2 = 64 2 resolved-scale Fourier modes, and 
N p = 512 2 subgrid-scale Gabor modes, and no additional viscosity. Run 
were also performed using a HDNS or the APVM model over 64 2 Fourier 
modes. The initial condition is a vorticity field with an energy spectrum 
concentrated at the forced wavenumber (k=15) and with a small amount 
of total energy. The simulation was forced by keeping constant the energy 
of the mode k=(15,0). In this situation, the vorticity field is progressively 
built via the stochastic forcing, which is itself strongly dependent on the 
exact structure of the vorticity field (since it must aim at keeping one mode 
constant). Due to the developing spectral cascades, the scale interactions 
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become less nonlocal and we may then expect sub-dominant terms to play 
an enhanced role (with respect to the decaying case). Indeed, we have ob- 
served that the vorticity fields in the DNS and in the model do not exactly 
correspond: both simulations display similar small scale intense vortices, 
but they are not located at the same places. This effect can also be seen 
more clearly in the spectra. They are shown in fig. |9[ Due to the chaoticity, 
they all differ at the mode k = 1 which is the most sensitive to the exact 
position of each individual vortex. At smaller scales, marked differences ap- 
pear between the models. Clearly, the HDNS gives the worst result, with a 
large deficit of energy over all scales. This is because in this case, the forcing 
scale is very close to the scale at which the dissipation takes place, and most 
of the energy is dissipated before the inverse cascade to larger scales can 
occur. The two other models, which do not introduce an explicit dissipa- 
tion at the cut-off scale, perform much better. The APVM model tends to 
overestimate the rate of energy at the largest scale and to underestimate 
it in the inertial range. Our model slightly underestimates the amount of 
energy at the largest scale, but gives results very close to the DNS in the 
inertial range of scale. At scales smaller than the cut-off, both the DNS 
and our model produce an energy spectra steeper than the fc -3 law. At the 
largest wavenumber, the beginning of the viscous range is clearly visible in 
the DNS but not in our model. This is again an indication of the large 
Reynolds number achieved by our model. As for computational times, the 
performances are similar to the decaying case. 
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Figure 8: Large-scale vorticity field (k < 21) after 50 turnover times in 
decaying turbulence, upper left: DNS on a 1024 2 grid and v = 1.8 1CT 4 , 
upper right: our model (M2 in table 2) with a separation scale at k = 21 
and 16384 modes in subgrid scales, lower left: simulation with hyper- 
viscosity (using u p k p io as the dissipation term with p=8) on a 64 2 grid and 
lower right: simulation with the APVM model (|23||) on a 64 2 grid. 
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Figure 9: Comparison of the energy spectra of forced turbulence after 150 
turnover times computed by different methods (our model refers to the run 
MO from the table 2) 
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5 Discussion 



We developed a new dynamic model of subgrid-scale turbulence based on a 
simple hypothesis about the subgrid-scale evolution. This approach is differ- 
ent from traditional turbulent models since our model provides expressions 
of the turbulent Reynolds sub-grid stresses via estimates of the sub-grid 
velocities rather than velocities correlations. The subgrid-scale dynamic is 
given by a linear equation, describing the advection of subgrid-scale wave- 
packets with the mean flow. This feature allows a reduction of the time 
step used in the simulation, via the use of a pseudo-Lagrangian method. 
We thereby achieved a reduction of the computational time by a factor 
150 in the typical cases we considered. Our method can also be used for 
Large Eddy Simulation strategies, in which only the large scales are com- 
puted. This allows a large memory savings, and an additional reduction of 
the computational time (typically a factor 10 in the case we considered). 
The resulting simulation is more costly than a traditional LES simulation, 
based on Fast Fourier Transform algorithms (like hyperviscosity or vorticity 
dissipative schemes). It could however become more competitive for more 
complicated geometries, where finite difference schemes become more ap- 
propriate than FFT. This would make LES models based on hyperviscous 
schemes much more costly, while our method would proportionally keep at 
the same computational performance. 

In the present paper, only 2D turbulence has been considered. As is well 
known, 2D turbulence is very special, because there is no vortex stretch- 
ing. One may therefore wonder whether our model could be applied to 



3D flows. To answer this question, Laval and Dubrulle 12] have recently 
performed a numerical analysis of the hypotheses pertaining our turbulent 
model. The main hypothesis is that the subgrid dynamics can be very well 
approximated by an equation in which non-local terms are retained, while 
local terms (non-linear in the subgrid velocities) are replaced by a turbulent 
viscosity. The numerical analysis, performed via a numerical simulation of 
the coupled resolved and approximate subgrid equations indeed showed that 
such hypothesis is valid, insofar as both the energy spectra, the structures 
and the statistical properties of such a model were very close to that of 
a DNS. Even so this analysis can only be performed at a rather moderate 
Reynolds number (80 to 200, based on the Taylor scale), we find this analysis 
encouraging. 
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A Derivation of the subgrid scale equation in (x,k) 
space 



Some properties of the GT will be useful in the sequel. They are (see [13, [T^] 
for details): 

(%u'(x,k,f) = «9;u'(x,k,i) (35) 
= ijfcju'(x,k,t) + 0(e*) (36) 

Uw'(x, k, t) ~ U(x, t)u'(x, k,t) + i (V x ■ V fc ) Uw'. (37) 
Let us derive the GT of the equation: 

d t uj' + div(UjOj') = F(x, t) + i/ t Ao/. (38) 

Using (|3"6|), we find that the GT of the viscous term is —vtk 2 Lu'. The GT 
and the time derivative commute, so that the GT of the first term of the lhs 
gives : 

d t u' = d t u>' (39) 

Using the space derivative property (p6|), the GT of the second terms can 
be developed into: 



djUju' = Ujdjuj' (40) 

d_ 

dki 

Ujdju' + idiUj— (ikjuty 



- UjdjW + idtUj—d^i 



dh 

d « 

~ Uj dj ijj 1 + i di Uj kj -^ru' — u'dj Uj 

OKI 

Using the incompressibility (djUj = 0), we finally obtain : 

D t u>(x, k, t) = F(x, t) - u t k 2 u)'(x, k, t) (41) 

with : 

D t = $ + U-V-V a! (U-k).V fc , (42) 

(43) 
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B Reconstruction formulae for the correlations 



Using the formula |36| for the GT of space derivative, one can derive the 
subgrid scale velocity field in the (x,k) space with respect to the subgrid 
scale vorticity: 

u'(x,y,p,q,t) = d x v'(x,y,p,q,t) - d y u'(x,y,p,q,t) (44) 
= d x v'{x,y,p,q,t) - d y u'(x,y,p,q,t) 
= ipv'(x,y,p,q,t) - iqu'(x,y,p,q,t) +0(e*), 

where k = (p,q) and x = (x,y). Reversing the formula, we obtain the 
expression of the Gabor Transform of the velocity: 



u'(x,y,p,q,t) = — — - uj'(x,y,p,q,t) + 0(e*), 
pz _|_ qz 

— ^p 

v'(x,y,p,q,t) = -5- — g v'(x,y,p,q,t) + 0(e*). 



(45) 
(46) 



Consider now the following expression: 
1 r 
2 



dk 



u'(x, k, t)u/(x, -k, t) + u'(x, -k, i)u/(x, k, t) 
u'(x, k, i)u/(x, — k, i) dk 
= J^t J /( e *(x - x'))e ik ( x - x V(x', t)dx' / /(e*(x - x"))e* k(x ~ x 'V(x", i)tfat" 

= y /(e*(x - x'))/(e*(x - x")K(x', t)w'(x", t) ^ e^'-^dk^ dx'dx" 

Using the definition of the Dirac function, 

1 /" e ^(x"-x') dk = 5(x /_ x '/ )) 
nY J 



dk 

(47) 



(2vr) 5 

and the fact that f 2 = G, one simply gets 
1 



u'(x, k, t)uj'(x, — k, i) dk = u'u'(x, t). 



(48) 



(49) 



We may proceed along the same line for the average of the non- linear product 
of a large scale field with a subgrid scale field. Using the definition of the 



average Uu/(x, t) can be written as follows, 

iy(x,t) = J / 2 (e^(x-x , ))U(x',ty(x , ,t)dx / . 



(50) 



29 



Using the Taylor development of U with respect to x at the first order, we 
obtain the first order approximation of this term: 

T3c7(x,i) = (UZ7)(x) +y'/ 2 (e ,t (x-x / ))(x-x , )VU(x / ,ty(x / ,t) ( ix / 
= (US7)(x) + 0(0. (51) 

If we now apply the definition of the average for a product of two subgrid 
scales field (|49|) with the quantities 1 et to', we finally obtain: 



Uw'(x,t) = U 
= U 



1 



(27T) 2 
1 



u/(x,k,i)i(x,-k,i) dk + (3(e*) (52) 
u/(x, k, i)/(A:)l rfk + 0(e*). (53) 



C Numerical computation of the Reynolds Stress 
components 

The subgrid scale field in physical space can be obtained from their discrete 
formula in the (x,k) by an integration with respect to k: 



1 



(2vr)2 /(0) 



u/(x, k, t)dk 



(54) 



2iV p 

X X Q (t)) <5(k-k Q (i))dk 



(2vr)2 /(0) 

y^y E £*( x ~ x "(*)) 
1 



/(0) 



X! ^a+(*)5'x(x-x Q , + (i))+ a- ai (t) 5* x (x - x a _(t)) > , 



a_=l 



where 



E=E + E> 

a=l a_i_=l a_=l 



(55) 



Sa+=i means sum over half the particles with the wavenumber k Q and 
J2a-=i ^ s * ne sum over particles with an opposite wavenumber — k Q . We 
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choose the "positive" particles a + with p a+ > and k a = (p a+ ,q a+ )- Be- 
cause x a+ = x a _, u/(x, i) can be written: 

w ' (x ' l) = 7(o) ^ (t) + * Q - 5x(x " Xq+ (56) 

Since u/(x, i) is real, the following property holds: 

*o+(f) = (57) 
and the formula used to rebuild the vorticity field in physical space is : 



2 ^ 

u/(x,i) = — J] 5x(x-x a+ (t)) 



(58) 



where 3? [o" a+ (t)] is the real part of <7 Q+ (i). The same developments can be 
made for the two velocity component, using ( f45|) and 



r3? [£<*+(*)] Sx(x-xa + (t)) (59) 



2 ^ -o 
u x,r = — r > — - — 

v ; /(0) p a+ 2 + 9 a+ 
^(x,*) = y^y E Pa /l + qa+ ^ ^W ^(x-x Q+ (t)) (60) 

where 3 s [<5" Q+ (t)] is now the Imaginary part of a a+ {t). 

The previous expressions can be used to compute velocity correlations. 
Using the analytical definition of these terms with respect the the subgrid 
scales field in (x,k) space (49,^), and using the definition of the subgrid 
scale discretization flITj|), we get: 



u'v'(x, t) 



1 



(2vr) 2 
1 



u'(x, k, t)v'(x, — k, t) 



dk 



(61) 



"IN* 



V %q a a 5 x (x - x a )5(k - k Q ) 



^ +ip 

K h p2+q2 



S x (x - x /3 )(5(-k - kp) 



dk. 
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If all the wave-packets have a different wavenumber (at least on a domain 
equal to the "support" of S x ), eq ^ can be written: 



N„ 



uV(x,i) = 2 R 

cr_l_ = l 



7<T a+ (T a _ S' x (x X, 



_0a+ 2 + <?a+ 2 ) 2 

Using the fact that a a _ and <r Q+ are complex conjugate, we finally get: 



(62) 



MM) = 2 £ 



<5"a4- 1 ^"x (" Xq,,), 



^ (.Pa+ 2 4" Qa + 



2 x 2 ra + 



/2 



(x,i) = 2 £ 



2^2 



(63) 



! (x,t)=2 y: 



(Pa + 2 + q a+ 2 ) 2 

The non-linear term involving the resolved scale velocity field can also be 



written in terms of wave-packets coordinates. Using eq. (53) and the defi- 
nition of discretization (|l6|) , one can write: 



Uw'(x,t) 



U 



(2vr) 2 
U 

u 

(2^F 



>(x,k,t)i(x,-k,t) 



dk 



u/(x,k,t)/(-k) 



dk 



2iV„ 



£ cr a /*(k)5 x (x - x a )£(k - k c 



a=l 



dk 



2U £ K [<7« + /*(k a+ )] £ x (x 

a_l_=l 



(64) 



D Choice of the filter 



The choice of the interpolating function S x (x) (eq. 17) dictates the shape 
of the filter and, therefore, of the unknown function f(x) in the definition 
of the Gabor Transform. The proof proceeds via the quantity u/ 2 (x, t), by 
equating the formula of the filter discretized on a regular grid with a cell 



32 



size Ax and the formula p8| 

+ \ _ sr,,> 2 



x,i) = ^^(xO/^x-x^lAx) 2 (65) 



= 2 |<r a+ | 2 S 2 (x-x a+ ) (66) 

o+=l 

For this equality to be valid at all grid point Aj, the following link between 
/ and S must hold: 

/ 2 (x) = C x <? 2 (x). (67) 
Using now the normalization: 

S'x(O) = 1. (68) 
J / 2 (x')dx' = 1. (69) 



the eq. |67j becomes: 
with 



/(x) = /(0)S x (x) (70) 



/(")- „ 1 (71) 

2/ ,/u ' 



5 x (x )dx 
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